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ABSTRACT 



In this paper we propose that additive self helicity, introduced by 



Longcope fc Malanushenko! (120081 ). plays a role in the kin k instability for com 



plex equilibria, simila r to twist helicity for thin flux tubes (iHood fc Priestl 11979 



Berger &: Field! Il984j ). We support this hypothesis by a calculation of additiv e 



self helicity of a twisted flux tube from the simulation of iFan fc Gibson! (120031 ) . 
As more twist gets introduced, the additive self helicity increases, and the kink 
instability of the tube coincides with the drop of additive self helicity, after the 
latter reaches the value of Ha/§ 2 ~ 1-5 (where $ is the flux of the tube and Ha 
is additive self helicity). 

We compare additive self helicity to twist for a thin sub-por tion of the tube 
to illu strate that Ha/& 2 is equal to the twist number, studied by lBerger fc Field 
(1l984l ). when the thin flux tube approximation is applicable. We suggest, that 
the quantity Ha/$ 2 could be treated as a generalization of a twist number, when 
thin flux tube approximation is not applicable. A threshold on a generalized twist 
number might prove extremely useful studying complex equilibria, just as twist 
number itself has proven useful studying idealized thin flux tubes. We explicitly 
describe a numerical method for calculating additive self helicity, which includes 
an algorithm for identifying a domain occupied by a flux bundle and a method of 
calculating potential magnetic field confined to this domain. We also describe a 
numerical method to calculate twist of a thin flux tube, using a frame parallelly 
transported along the axis of the tube. 
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Introduction 



According to a prevalent model coronal mass ejections (CMEs) are triggered b y current- 
driven magnetohydrodynamic (MHD) instabili ty related to the external kink mode (IHood &: Priest 
19791 : iTorok et allhooi iRachmeler et alil2009h . The external kink mode, in its strictest form, 
is a helical deformation of an initially symmetric, cylindrical equilibrium, consisting of heli- 
cally twisted field lines. The equilibrium is unstable to this instability if its field lines twist 



abou t the axis b y more than a critical angle, typically close to 3n radians (IHood fc Priest 



19791 ; iBatyl l200ll ) . The helical deformation leads to an overall decrease in magnetic energy, 
since it shortens many field lines even as it lengthens the axis. 

Equilibria without symmetry can undergo an analogou s form of current-driven instabil- 
ity u nder which global motion lowers the magnetic energy (IBernstein et al.lll958l ; iNewcomb 
19601 ). Such an instability implies the existence of another equilibrium with lower magnetic 
energy. The spontaneous motion tends to deform the unstable field into a state resembling 
the lower energy equilibrium. Indeed, it is generally expected that there is at least one min- 
imum energy state from which deformation cannot lower the the magnetic energy without 
breaking magnetic field lines; its energy is the absolute minimum under ideal motion. 

Linear stability and instability are determined by the energy change under infinitesimal 
motions. An equilibrium will change energy only at the second order since first order changes 
vanish as a requirement for force balance. Ideal stability demands that no deformation 
decrease the energy at second order, while instability will result if even one energy- decreasing 
motion is possible. The infinite variety of possible motions make it impractical to establish 
stability in any but the simplest and most symmetric equilibria. 

Based on analogy to axisymmetric systems it is expected that general equilibria, includ- 
ing those relevant to CMEs, are probably unstable when some portion of their field lines are 
twisted abou t one another by more than some critical angle. This expectation was mentioned 
in a study by lFan fc Gibson! (120031 ) of the evolution of a toroidal flux rope into a pre-existing 
coronal arcade. They solved time-dependent equations of MHD in a three-dimensional, 
rectangular domain. Flux tube emergence was simulated by kinematically introducing an 
isolated toroidal field through the lower boundary. The toroidal field was introduced be- 
neath a pre-existing arcade slowly enough that the coronal response never approached the 
local Alfven speed. Fan and Gibson concluded that the system underwent a current- driven 
instability after a critical amount of the torus had been introduced. They bolstered this 
claim by performing an auxiliary run where the kinematic emergence was halted and the 
system allowed to evolve freely; it settled into an equilibrium. 



While twist angle has proven useful in a few cases, it is difficult to demonstrate its utility 
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as a threshold in general, asymmetric equilibria. Indeed, in any but a few very symmetric 
cases there is no simple, obvious way to define the angle by which the field lines wrap about 
one another. The local rate of twist is given by the current density, which is after all the 
source of free energy powering the instability. On the other hand, excessive local current 
density is not sufficient to drive instability. This fact is illustrated by numerous examples of 
discontinuous field which are minimum energy states. 

It has been suggested that a threshold exist s, in general equil i bria, for s ome global 
quantity such as free magnetic energy or helicity (jZhang et al.ll2006l ; iLowi Il994j ) . If this is 
the case then we expect the instability to lower the value of this global quantity so that it 
falls below the threshold value in the lower-energy, stable equilibrium. Magnetic helicity is 
a logical candidate to play this role since it is proportional to total twist angle in cylindrical 
fields. Relative helicity in particular is a proxy for currents. Helicity is, however, conserved 
under ideal motion and therefore will not be reduced to a sub-threshold value by an ideal 
instability. 

The total helicity of a thin, isolated flux tube can be written as a sum of two terms 
called twist and writhe faerger fc Field! ll984J : IMoffatt fc Riccalll992h . 



H = Ht + H\y. 

The writhe depends on the configuration of the tube's axis while the twist depends on the 
wrapping of field lines about one another. A cylindrical tube has a perfectly straight axis and 
therefore zero writhe helicity. Any ideal motion which helically deforms the entire flux tube 
will increase the magnitude of the writhe helicity. Since the motion preserves total helicity 
the change in writhe must be accompanied by an offsetting change in twist helicity. If the 
writhe has the same sign as the initial twist, then the motion will decrease the twist helicity. 
In cases where th e magnetic energy depend s mostly on twist, this motion will decrease the 
magnetic energy (ILinton &: Antiochosl 120021 ) . The straight equilibrium is therefore unstable 
to an external kink mode. 

Topologically, the foregoing properties of magnetic field lines could be compared to the 
properties of thin closed ribbons. One may introduce twist number, writhe number and their 
combination, called linkage number, is a pre served quantity in the absence of reconnection 
faerger & Field! Il984t IMoffatt & Riccalll992h . 



L = Tw + Wr. 



By analogy to the case of a thin isolated flux tube we consider the twist helicity, rather 
than the total helicity, to be the most likely candidate for a stability threshold. Indeed, 
within a thin flux tube it is possible to derive a net twist angle among field lines and 



-4- 



H T = $ 2 Tw = $ 2 A6»/2vr, where $ is the total magnetic flux through a cross-section of the 
tube and A8 is the net twist angle. 

Twist and writhe are, however, defined only in cases of thin, isolated magnetic flux 
tubes, and can no more set the threshold we seek than the net twist angle can. 



Recently iLongcope fc Malanushenkol (120081 ) introduced two generalizations of relative 



helicity applicable to arbitrary sub- volumes of a magnetic field. They termed both general- 
ized self-helicity, and the two differed only by the reference field used in their computation. 
The one called additive self-helicity (that we denote Ha) uses a reference field confined to 
the same sub-volume as the original field, and can be interpreted as a generalization of the 
twist helicity to arbitrary magnetic fields. The additive self-helicity of a thin, isolated flux 
tube is exactly the twist helicity. 

Since the additive self-helicity can be computed for arbitrary magnetic fields we propose 
that it (normalized by the squared flux) is the quantity to which current-driven instability 
sets an upper limit, which could be considered a generalized twist number: 

Tw (gen) = H A /&. (1) 

The paper is organized as follows. In Section 2, we describe a method for calculating 
additive self helicity and Tw^ gen ) numerically. There are two large and nontrivial parts of 
this calculation, that we describe in 2.1 and 2.2: locating a domain containing a given flux 
bundle and constructing a potential field in this domain by Jacobi relaxation. In Section 
3, we apply t he method to a simul ation to support our hypothesis, the emerging twisted 



flux tube from lFan &: Gibson! (120031 ). In 3.1 we briefly describe this simulation, and then in 



3.2 we show different embedded do mains defined by differe nt subportions of the footpoints. 



In 3.2 we describe, how the twist of iBerger fc Field! (119841 ) could be calculated for those of 



the domains for which thin flux tube approximation is applicable. In Section 4 we present 
the evolution of additive self helicity, unconfined self-helicity, twist (for "thin" domains) 
and the integrated helicity flux in the simulation. We demonstrate that Tw( gen ) increases 
corresponding to helicity flux, that it drops after it reaches a certain value (about 1.5) and 
that this drop coincides with the rapid expansion of the tube due to the kink instability. We 
also demonstrate that the unconfined self helicity grows only when helicity flux is nonzero and 
that it stays constant when kink instability happens. We also show that Twi gen \ corresponds 
to Tw when thin flux tube approximation is applicable. 
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Numerical Solutions 



The object of study is a magnetic field B (r) defined in a domain P, r G D, that lies on 
and above the photosphere, z > 0. By domain we understand a volume that encloses the field: 
B • n = on all boundaries, &D, except at the photosphere, where B • n = B z (x, y, z — 0). 
An example of such a volume is the coronal part of an ^-shaped loop. The self-helicity is 
given by 

H A (B,P(V),V) = J (B-P)-(A + A P )dV, (2) 
v 



as defined in lLongcope &: Malanushenkol (120081 ) . Here P is the potential magnetic field, 
whose normal component matches the normal component of B on the boundary dT>, 



■ n 



B ■ n 



&D, 



(3) 



A and A P are the vector potentials of B and P respectively (as discussed in Finn fc Antonsen 
(119851 ). helicity, defined this way is gauge- independent). 



Once the self-helicity is known, the twist is given by eq. ([I]) with $ being the total 
signed flux of the footpoints of the configuration: 



$ = J B z dxdy = - J B z dxdy. (4) 

z=0,B z >0 z=0,B z <0 

In the next two sections we discuss methods of numerically obtaining V, from given foot- 
points, and P. 



2.1. Finding the domain. 

In order to describe the domain on a grid we introduce the support function: 



1, if r G V 
0, if r i V. 



This is a function of the given magnetic field B and some photospheric area, called the 
boundary mask. By definition, every field line, initiated at any point on the boundary 
mask and having the other footpoint somewhere within the mask, is completely inside the 
domain T>. If the field line traced in both directions from some coronal point ends within the 
photospheric mask, then this point also belongs to the domain. In numerical computations 
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we replace "point" with a small finite volume, voxel Vijk (3-dimensional pixel). We define a 
voxel to be inside T> (equivalent to saying O (i*^) = 1), if there is at least one point inside 
it that belongs to V. 




Fig. 1. — An example of what we call a domain. Here the field from four photospheric sources is 
computed on a half-space. Two possible domains are shown in two different colours. 

The simplest method of constructing the support function would be to trace a field line 
in both direction from every voxel of the computational grid, set G = 1 in the voxel if the 
footpoints both terminate in pixels from the boundary regions, and set = otherwise. 
This, however, is a very time-consuming algorithm, especially for a large arrays of data. In- 
stead we use an algorithm which reduces the computational time by tracing field lines from 
a subset of voxels. It works by progressively adding voxels to adjacent to those already 
known to belong to T>. 

We add a voxel centered at r it j^ to the domain under two different circumstances. 1. 
A field line initialized somewhere within the volume of the voxel Uij,jfe, centered at r^j^, is 
found to have both footpoints within the boundary mask. 2. A field line initiated in some 
other voxel, and determined to belong to £>, passes through some portion of the volume v^j^. 

Initially, the domain consists only of footpoint voxels, so the initial step is to trace field 
lines initiated at the footpoints, assuming, that at least some of these lines will lie in the 
domain. 

We illustrate the method on a simplistic case of a potential magnetic field, confined to 
a half-space, with B z = everywhere at the photosphere, except at four pixels, as shown 
on Fig. [2j We have computed the magnetic field inside a small box of 15 x 15 x 15 pixels, 
centered around the photospheric sources. The boundary mask consists of these four voxels 
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at the photosphere with non-zero vertical magnetic field. In this simplistic example the 
initial guess would be four field lines, initiated at four footpoint voxels, as shown on Fig. (2J 
left (note that in this particular example a field line, initiated at one voxel, ends at another 
voxel within the mask and thus is the same as the field line, initiated at that another voxel, 
so these four initial guesses are really two, not four field lines). The voxels of the initial guess 
are shown with crosses. 

In an algorithm, this would be the first step: 

Step 1: Make the initial guess: trace field lines from the footpoints. 

As soon as an initial step is made, the next step is to assume, that the immediate neigh- 
bourhood of voxels known to be in T> are likely to be also in the domain. Thus, in the next 
(iterative) search the following steps are performed: 

Step 2: Locate voxels on the boundary of the current domain. 

Step 3: For every voxel on the boundary: trace a field line and check whether it is in the 
domain. 

If yes: Add the voxel to the domain. Add all voxels along the line to the domain. Exclude 
them from the boundary (there is no need to check them again). 

If no: Mark the voxel as "questionable" . (If there is a field line, which passes through 
the voxel and does not belong to the domain, then at least part of the voxel is 
outside of the domain. Since its immediate neighbourhood is in the domain, then 
it is possible that part of it is also in the domain.) 

Loop: Repeat steps 2-3 until all the voxels in the boundary are "questionable" and no new 
voxels are added. 

When the iterative search does not find any new voxels, we make the final check of 
the boundary voxels. The idea is to trace field lines from all corners of such "questionable" 
voxels to see, which corners (and thus which part of a voxel) belongs to the domain. We 
consider this to be optional check, which may improve the precision of the definition of the 
domain by at most one layer of voxels. 

This last search may also give information about the normal to the domain surface. If 
it is known that some corners of a voxel are in the domain and some are not, it is possible 
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to approximate the boundary as a plane separating those two groups of corners. 



Step 4, optional: For each voxel, marked previously as "questionable", check the corners (by tracing field 
lines) to see which of them are in the domain and which aren't. Keep this information. 



-' 
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Fig. 2. — fie/tj — The first iteration of the iterative search: from the initially selected voxels 
(crosses), check those surrounding (circles) for membership in the domain. Repeat until no "sur- 
rounding" voxels can be added to the domain, (middle) — The voxels, checked on all iterations in 
the middle plane. For every field line, a cross shows where it was initialized. Yellow are "accepted" 
lines (and thus all voxels that contain them are "accepted") and black are "not accepted" lines 
(and thus only voxels where these lines were initializes from are "not accepted"), (right) - The end 
result. The green crosses mark voxels that are found to belong to the domain and the green circles 
are the neighbourhood of the domain. White (initial), yellow (iterative) and red (final) field lines 
are traced and found to be in the domain; black lines are found to be not in the domain. Note that 
the domain is "covered" by much fewer lines than an exhaustive search would do. 
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2.2. Constructing the Confined Potential Field P 

Once the domain has been determined, the next step is to construct the potential mag- 
netic field confined to it. We use a common relaxation method on a staggered grid in orded 
to account for the complex boundaries of V. 



We introduce a scalar potential B p = V% and look for the solution of the Laplace's 
equation for x 

V • Bp = V 2 x = 0. 

By the definition of T>, field lines never cross dT>, except at the lower boundary, z — 0. 
Thus, boundary conditions for B p could be written as: B p ■ n\ dv z ^ = and B p ■ n| 9x , z=Q = 
H z (x,y). This is equivalent to Neumann boundary conditions for \'- 



dx I 

dn I dV,z^0 
&X I 



0. 



dz\ov,z=o B z (x,y). 
The Algorithm for the Relaxation Method 



(5) 



We use the Jacobi iterative method (see, for example iLeVequel Il955l ) to solve for the 
potential field. Here we briefly summarize the algorithm and further explain in details. The 
n + 1-th iteration is 



1. Vr G D: calculate a new iteration ^[ n+1 l as a solution of the equation ^[ n+1 l — — 
Kh 2 \7 2 x^ n \ where h is the grid spacing. The Laplacian V 2 x'™' (r), found ausing stan- 
dard finite difference methods, is equivalent to an average over some stencil of neigh- 
bouring points minus the central value; if is a constant that depends on the exact 
shape of the stencil. 

2. Vr 6 G dV: set x [n+1] (r 6 ) so as to satisfy boundary conditions (BCs). 

3. Repeat steps 1-2, until the difference between x^ ( r ) and x^ n+1 ^ (r) is sufficiently small 
in some sense (namely, until ||x' n+1 ^ — X^W < e ; where e is pre-defined small number). 



Staggered Mesh 



The functions B x (x, y, z), B y (x, y, z) and B z (x, y, z) are defined on the same mesh points 
(xi, yj, z k ). If we are interested in finding x(x, y, z), so that B x = B y = |^ and B z = 
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it advantageous to define x ^ n between the original mesh points and calculate the derivatives 
using finite difference as following: 

, . _ Xfa+X/2, Vj, Zk) - Xfa-l/2, Vj, Zk) 

^x\ x i,Vj, z k) — > J 

x i+l/2 ~ x i-l/2 

and so on for B y and B z . x, then, would only be defined in the middle of the faces of cubic 
voxels, i.e., at points (i ± 1/2, j, k), ± 1/2) and (i,j,k± 1/2). 

Such a mesh, called a "cartesian staggered mesh", is known to have better numerical proper- 
ties, such as immunity from decoupling of variables and having a smaller numeric dispersion 



(jPerotll2000l . see, for example). 



The finite difference approximation of a Laplacian at one point can be interpreted as a 
weighted average over a stencil of several points minus the value at that point. For example, 
in the 2D case the second order approximation to V 2 x{%, y) on a uniform Cartesian grid at 
the point fa, yj) could be computed over a 5-point stencil: 

V 2 x (xi, w — (x (xi-i, yj) + x fa+i, Vj) + X fa, Vj-i) + X fa, Vj+i) - 4% fa, y d )) 

(here h is the spacing of the grid). It could be rewritten as 

1 h 2 

X fa, yj) ~ - (x fa-i, Vj) + X fa+i, Vj) + X fa, Vj-i) + X fa, Vj+i)) - -J^ 2 X fa, Vj) ■ 

The Jacobi method uses this equation to iteratively update the value at the point, constantly 
assuming V 2 x = 0. In the case of the 5-points stencil the updated value would be 

X [n+1] fa, Vj ) = \ {x [n] fa-x, Vj) + X [n] fa+i, Vj) + X [n] fa, Vj-i) + x [n] fa, y j+ i)) ■ 

In our case of a 3D staggered mesh, choosing a stencil becomes more complicated. We 
propose a 13-point scheme, shown on the right of Fig. [3] (black dots). To motivate this stencil, 
we derive it from the "unstaggered" one (Fig. [3], left, gray dots). In an "unstaggered" finite 
differencing scheme the [n + l]-th iteration in Jacobi method would be expressed as 

Qx [n+1] (0) = X [n] (A) + X [n] (A 2 ) + x N (5i) + X [n \B 2 ) + X [n \Ci) + X [n] (C 2 ). 

But for the staggered mesh x is undefined at these nodes. This can be resolved by setting x 
at each "gray" point to be equal to the average of its 4 closest neighbours, 

X^(Ai) = \ [x^KSA,) + x [nl (7Mi) + X [n] (SA 1 ) + X [n] (0)} , 
X [n] (Bi) = \ [x [n] (SB 1 ) + x H (TB 1 ) + x [n KSBt) + x [n] {0)] , 
X [n] (Ci) = \ [x [nl (^0 + X H (TA 2 ) + x^KTBj + x [n] (TB 2 )] 
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and so on. Then we may substitute this in the original expression and get: 

6x [n+1] (0) = 2 x \ [x [n] ( T Ai) + X [n] {TA 2 ) + x ln] { B ^i) + x [n] {BA 2 )] + 
+ 2 x i [x [n] {TB x ) + x [n] {TB 2 ) + x H (BBi) + X H {BB 2 )] + 
+ | [ X [n] {SA x ) + x [n \SA 2 ) + X {n] {SB x ) + X W (^ 2 )] + 
+ 4x± x M(0), 

which is eqivalent to 

X [n+1] (0) = i [x W (^i) + X W (TA 2 ) + x W (5^i) + X W (^ 2 )] + 

+ i [x w (^0 + x w (^ 2 ) + x w (55i) + x w (55 2 )] + 

+ £ [x W (^i) + X N (5A 2 ) + x W (55i) + x w (^ 2 )] + 

+ k w (0). 

With these weights th "farthest" nodes 5[AS]r 12 ] have half the influence on the laplacian, of 
the "closer" nodes. Note also, that the sum of the weights is one. 




Fig. 3. — The averaging kernel for the laplace's equation on 3D staggered mesh (right) and the 
motivation for it (left). For example, the stencil for a face with normal vector z would include five 
"z faces" (including itself), four "x faces" and four "y faces" (two of each above and below). 

Boundary Conditions 



Boundary conditions (given by Eqn. [5]) in the staggered mesh is particularly easy if 
one assumes that the boundary surface passes inside of boundary voxels, rather than on 
their sides. Suppose, for example, that the boundary plane normal to z passes through the 
center of the voxel v^. Then the BC for this voxel would be that B z k) = 0, or simply 
X (hj,k+ \) =X {i,j,k- \). 
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To motivate such choice of the boundary, we note that boundary voxels, by definition, 
are the voxels part of which is inside of T> while part is outside. Such a conclusion is made 
about voxels, some of whose corners are inside of T>, and some of the corners are outside of 
T> (this information about the domain is obtained in the step 4 of the algorithm, described 
in section [2~T|) . We approximate the boundary inside of each boundary voxel as a plane, 
that passes through the center of the voxel and that separates its "exterior" part from its 
"interior" part. Such approximation will err by no more that 1/V2 voxel's length off the 
real location of the boundary. We also find it easier to work in terms of faces rather than 
corners, since this is where % is defined. (We say, that a face is "exterior" to the domain if 
more than two of its corners are not in the domain, i.e., for a voxel, we say, that if only one 
corner or only one edge are "exterior", we do not consider it a subject to BC's). 

There are several ways to orient such a boundary plane inside a voxel, based on the 
behaviour of the boundary in the immediate surrounding of the voxel. 

1. The voxel has only one face outside of the domain. Then we consider the boundary 
parallel to that face of the voxel (see Fig. HJ left). If, say, the boundary is parallel to 
the face between faces A and A\ (see Fig. HJ bottom left), then the normal field to the 
boundary is B • AA X (hereafter AA X denotes a unit vector along the line from A to 
Ax, which might be ±x, ±y or ±z), and BC would be formulated as 

Xa = l-XA' + O-Xb' + • xc ■ 

2. The voxel has two adjacent faces outside of the domain. Then we approximate the 
boundary as a plane, that cuts off these two faces, as shown on Fig. HI middle. If faces 
A and B are outside and faces Ax and B\ are inside of the domain, then we consider 
the normal field to be B • 4= ( AAi + BBi ) and set BC's as 



Xa = ■ xa' + 1 • Xb> + ■ xc, 
Xb = 1 • Xa> + O-Xb' + ■ xc- 

3. Similarly, if three mutually adjacent faces of the voxel are outside of the domain (and 
three others are inside), as shown on Fig. HI right, then, analogously, we assume that 
the normal field is B • ^ ^AA X + BBx + CCi j and BC's could be set in the following 

way: 

Xa = 0-XA' + \-Xb> + \-Xc, 
Xb = \-Xa< + 0-xb' + I'Xc, 
Xc = \-Xa> + \-Xb> + O'Xc- 

(Note that in this case there are really three variables and one equation to satisfy; 
thus, there are different solutions to x- But each of those solutions would be valid, as 




long as it satisfies B • n = 0.) 



4. "Everything else": the voxel has three or more non-adjacent faces that are outside of 
the domain, but still is on the boundary. It is considered an extraneous voxel and is 
removed from the boundary. 




Fig. 4. — Different ways to approximate the boundary surface inside of a boundary voxel, depending 
on which portion of the voxel is found to be outside of the domain. White dots are the centers of 
the "interior" faces, gray dots are the centers of the "exterior" faces (see explanation in the text), 
the thick plane is the proposed approximation of the boundary surface dV. 
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3. The Experiment 

The method described above was tested on a simple quadrupole exam ple, and the values 



of sel f-helicity it gives are in a good agreement with theoretical predictions (jLongcope fc Malanushenko 



20081 ). That work, however, does not consider any sort of stable equilibrium and does not 



study any kinking instability thresholds, similar to those developed in Hood & Priest, 1981. 



The objective of the current work is to test whether the parameter Ha/Q 2 behaves like 
a total twist in the sense that it has a critical value above which a system is unstable to 
a global disruption. To do so, we use the numerical simulation of kink instability in an 
emerging flux tube from Fan & Gibson, 2003 . 



3.1. Simulation Data 



The initial configuration is a linear arcade above the photosphere, into which a thick, 
non-force-free torus was emerged. Inside the torus the field lines wind around its minor axis 
and the field magnitude drops with distance from the minor axis. The exact shape of the 
magnetic field, in the coordinates shown on Fig. [51 is the following: 

B = + B v <p = B t e~™ 2/a2 (q^# + , (6) 

where a = 0.1L is the minor radius, R = 0.375L is the major radius, q = —1, B t = 9B , L 
is the length scale of the domain (further in our calculations L — 1), Bq is the characteristic 
strength of the photospheric arcade the torus is emerging into, and the time is given in the 
units of Alfven time, ta = L/va- The field strength drops as e~ ro2//a2 with w being the 
distance from the minor axis. At w = 3a magnetic field was artificially set to 0. 



The torus is "emerged" from underneath the photosphere with a constant speed. There 
is a mass flow across the photosphere in the area, and the emerging tube is driven into the 
domain by an electric field at the boundary. This is made in the following way: for each 
time step (starting at t=0 and until the axis of the torus has emerged, t=54) the vertical 
photospheric field is set to that from the appropriate slice of the torus's field. Dynamical 
equations are then solved in order for the field above z = to relax, so that at every time 
step the resulting configuration is a force-free equilibrium. The unsigned photospheric flux 
clS 8u function of time is shown on Fig. 
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Fig. 5. — B is set in spherical coordinates (r,9,(j)) with the polar axis directed along y. We will 
mainly use different coordinates, namely, (zo,ip,(f>). R is the major radius of the torus, a is the 
minor one, p = rsin(#) is the distance from the y axis. 

A visual representation of characteristics times is shown on Fig. Different rows cor- 
respond to different times: t = 15 - the tube is about to emerge; t = 24 - the minor axis 
of the torus has emerged; t = 32 - the bottom of the torus has emerged; t — 45 - the tube 
undergoes acceleration; t — 54 - the major axis of the torus has emerged, the torus has 
stopped emerging, the tube starts getting a significant writhe; t — 58 - the tube escapes 
the domain; the simulation is over. Note that the torus starts to kink at t > 45 and keeps 
kinking until it escapes the computational domain at t = 58. 



*(t) 
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Fig. 6. — The total amount of unsigned photospheric flux, as defined in Eqn. ([6]), of the torus 
(not counting the arcade), plotted as a function of time. The major axis of the torus emerges at 
t = 54. The maximal value of flux is reached earlier than that because of the field winding aroung 
the torus and thus being not necessarily normal to z. After t = 54 the torus has stopped emerging 
and thus the magnetic field at the photosphere remains constant. 
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Fig. 7. — The characteristic times for the simulation of Fan & Gibson, 2003, different rows corre- 
spond to different time (see detailed explanation in text). First column - XZ slices, the analytical 
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3.2. Computing Ha For Given Volume And The Potential Field. 

We define different domains, T>, with the same field, by making a different choices of 
boundary mask. We were interested in how different portions of the torus, namely, the "core" 
and the outer layers behave during the instability. 

Our masks are defined to be within the photospheric intersection of the emerging torus, 
ro < ttmax- By choosing different values of w max we construct domains, containing different 
portions of the emerging flux tube The footpoints of domains with different w max are shown 
on Fig. [HJ The shape is distorted with respect to the original cross-section of a torus due to 
reconnection with the arcade, current sheet formation and due to near horizontality of some 
field lines. 

We found domains for masks with w max G [0.5, 1.0, 2.0] R at different times during the 
emergence. We computed 6 {w maxi t) and then constructed a potential field confined to it. 
The results are shown in Fig. [BJ Fig. [9] and Fig. [TUl 




Fig. 8. An exampie of footpoints of domains (w max 6 [0.5, 1.0, 2.0] , t = 50). The verticai field, 
B z , is shown in grayscale and horizontal field is shown with arrows. Three pairs of concentric curves, 
counting from inside out enclose footpoints of the domains defined by vj max = 0.5, w max = 1.0 and 

t^max — 2.0. 
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Fig. 9. — (left column) - field of a torus, confined to domains of different vj max , with footpoints 
shown on Fig. [51 (right column) - the potential field, constructed for each such domain. 




Fig. 10. — The strip plot of the results of the computation. The original data is shown above 
and the relaxed notential field P - below the nhotosnhere. The dotted line indicates slices of the 
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For each t and w max we calculated vector potentials of the actual field, 9 (t, w max ) B (r, t) , 
and the reference field P (r, t, w max ). To do this we used a gauge in which one of the com- 
ponents of the vector potential (in our case, A z ) is identically zero. The other two could be 
found with a straight-forward computation: 

z 

A x (x,y,z) = J B y (x,y,z')dz' 
o 

z 

Ay(x,y,z) = f(x,y) - J B x (x,y,z')dz' (7) 

o 

x 

f(x,y) = J B z (x',y,0)dx' 
o 

In terms of these elements the addirive self helicity the additive self-helicity: 

H A (t,w max )= j (&B-P)-(A + A P )dV (8) 

0(t,ro max ) 

is computed. 



3.3. Measuring Twist in Thin Flux Tube Approximation 

To make contact with previous work we compare the additive self helicity to the twist 
helicity in our flux bundeles. It can be shown analytically that in the limit of a vanishingly 
thin flux tube these quantities are identical. Here we must compute twist helicity for flux 
bundles of non-vanishing width. We do this in terms of a geometrical twist related to twist 
helicity. 

One cannot really speak of twist, or of an axis, in the domains defined above. First, 
the thickness and the curvature radius of the flux bundles are comparable to their lengths. 
Secondly, the magnetic field and the twist vary rapidly over the cross section of the bundle. 

The domains constructed from the smaller masks, w max = 0.5 and zu max = 1.0, may, 
however, be suitable for approximation as thin tubes. Even in these cases the approximation 
may suffer near the top part at later times: at t = 50 the radius of curvature becomes com- 
parable to the width, and later, during kinking the radius of the tube becomes comparable 
to the length (see Fig. [H]and Fig. [7]). 

We define an axis for the flux bundle by first tracing many field lines within it. Then 
we divide each field line into N equal segments (N is the same for all lines) of length Li/N, 
where L, is the length of the i th line. If the bundle were an ideal cylinder, the midpoints of 
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the n segment from every line would lie on a single plane; provided the bundle is thin the 
these midpoints will lie close to a plane. We define the n th point on an axis by the centroid 
of these approximately co-planar points. The set of N centroids forms the axis of our tube. 

We then define the tangent vector li along this axis, and a plane normal to this vector 
and thus normal to the flux tube (at least in the thin flux tube approximation). If the tube 
has some twist in it, then the point where one field line intersects the plane will spin about 
the axis as the plane moves along the tube. Such spinning must be defined relative to a ref- 
erence vectore on the plane which "does not spin" . The net angle by whcih the intersection 
point spins, relative to the non-spinning vecotr, is the total twist angle of the tube. In a thin 
tube all field lines will spin by the small angle; in our general case we compute an everage 
angle. 

We produce a non-spinning reference vector using an orthonormal triad, arbitrarily 
defined at one end of the tube, and carried along the axis by parallel transport. For a curve 
with tangent unit vector I, the parallel transport of a vector u means /• (du/dl) = 0. To 
impliment this numerically an arbitrary unit vector u is chosen at one end of the axis 
perpendicular to the tangent, uo ■ k> = 0. The third member of the triad is vo = uo x Iq. At 
the next point, U\ is chosen by projecting u onto a plane normal to l\ and normalizing it 



(see Fig. [TTT) . Then vi = Ui x Zi, and the procedure is repeated for every segment along the 
axis. 



Based on the analogy between Tw and Ha/ $ 2 , it would be natural to introduce quantity 
analogous to L and Wr in a similar way. We propose that L in the general (non- "thin" ) case 
might be analogous to the unconfined self-helicity, introduced in Longcope & Malanushenko, 
2007, and Wr is similar to the helicity of the confined potential field relative to the unconfined 
potential field. From equation (3) of Longcope & Malanushenko, 2007 




4. 



Results 
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Fig. 11. — An illustration of parallel transport of a coordinate system. At every next step one 
unit vector of the previous coordinate system is projected to a new normal plane and normalized; 
the second vector is created anew as perpendicular to the new unit vector. I is the tangent vector 
of the axis, n is the unit vector in normal plane, carried with the plane along the axis. 



X 

H{B, P v , V) = J d 3 xB ■ A - J d 3 xP v • A P + J dxdyB z (x, y, 0) J dx [A(x') - A P (x')] 

V V z=0 x 

(9) 

(where x = r(x,y,0) and Pv is a potential field confined to V that matches boundary 
conditions B • n\g V = P v • n|ay) by plugging it into H(Ox>B, P©, V) and H(B, Pv,e, V) and 
adding them together it immediately follows, that 

H (Q V B, Pv, V) + H (P v , P v ,e, V) = H (B V B, P v ,e, V) , (10) 

where T> C V and is a support function of T>. By P© we mean the potential field con- 
fined to T> (and identically zero outside of T>) that matches boundary conditions B - n|gx> = 
Px> ■ n|az3, and by Py,e we mean the potential field, confined to V, that matches bound- 
ary conditions Ox>B • n|gv = Pv,e ■ ^\dv- As long as T> is fully contained in V, which is 
constant in time, the quantity iJ unCi v/ < £> 2 = H (Ox>B, Py,e, V) /$ 2 will behave like L and 
H (Px>, P v ,e, V) /$ 2 would then behave like Wr. 



Fig. [TH compares the generalized twist number, if^/$ 2 , with helicity, unconfined to the 
flux bundle's volume, but confined to the computational domain of the simulation: -ff unC) box- 
In this case V is the computational domain, a rectangular box. The behaviour of all quantities 
matches expecation: H unc ^ ox /Q 2 increases as the torus emerges, and stays nearly constant 
after the emergence is complete (the slight decrease is due to the reconnection with the 
arcade field). The generalized twist number, Ha/Q 2 also increases with the emergence, but 
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Fig. 12. — An example of the ctxis, found, for TU max = 0.5, t = 58, and the corresponding coordinate 
system, carried along by parallel transport. I, u and v are drawn in red, green and blue colors 
respectively. 

decreases between t = 50 and t = 58 - the time when the torus kinks (see Fig. [7]). For 
different zu max the decrease seems to start at a slightly different time. 



Fig. fHl demonstrates as well, that the general behaviour of H mu ~y/$ 2 is qualitatively 
similar whether the volume V over which unconfined helicity is computed is the compu- 
tational domain or the half space. To compute the unconfined helicity in th e half space, 
H unr 7, v , we integrate the helicity flux in the way described in flDeVorel 120001 ) and used in 
(IFan &: Gibson! 120041 ). The helicity flux is computed relative to the potential field in half 
space, and thus, the helicity flux, obtained in this way, might be considered a "confined to 
a half space" . 

Longcope and Malanushenko (2008) show that i? U nc,box = H unc ^z + when the volumes, V 
and Z + and the vertical field, B z (z = 0), all share a reflectional symmetry. This situation 



Fig. 13. — An illustration of how kinking decreases twist. An axis (solid) of a "thin", w max = 0.5, 
tube and a single field line (dotted with diamonds) at a different times: top row is t = 50, the 
field line has A8 ~ — 3.17T and bottom row is t = 58 and the field line has A9 ~ — 2.4-7T (and 
Tw = A6/2ir). Left column is sideview and right column is the trajectory of the line in the tangent 
plane with coordinate system decribed above. 

occurs in the simulation only for t > 54 when the torus is fully emerged and its major axis is at 
the photosphere. At these times the vertical component of the field is the toroidal component 
of the torus, which is symmetric about y = 0. Due to reconnection with the arcade, however, 
the footpoints of T> may not share this symmetry, in which case the photospheric field QvB z 
is not precisely symmetric. If the two helicities were ever to coincide, it would be at t = 54, 
so we choose constant of integration by setting i? unc ,box = H unc ^ z+ at that time. The time 
histories of both unconfined helicities are plotted in Fig. [TU The discrepancy between the 
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two before t = 54 arises from the non-vanishing helicity of Py,e relative Pz + ,e owing to a 
photospheric field, B z , lacking reflectional symmetry. In spite of the discrepancy, we draw 
from each curve the same basic conclusion, that the kink deformation of T> does not change 

H un c,V- 



H/<t> 




30 35 40 45 50 55 

:ir"e 



Fig. 14.— The comparison between Ha confined to the volume of the flux tube), H unc \ )OX 
(confined to the box in which the original simulation was performed) and H unc z + (confined to half- 
space), normalized by $ 2 . Vertical dashed line at t = 54 indicates the time when the emergence 
has stopped and all further changes in Tw would be due to kinking and numerical diffusion, and 
all earlier changes are altered by the emergence of the tube and thus non-zero helicity flux over the 
surface. For w max of 2.0 and 1.0 it's clearly visible, that: a) after t = 54 the unconfined helicities 
remain nearly constant, while the confined to flux bundle's T>, that is, additive self helicity, decreases 
due to kinking; b) before t = 54 the difference between H unCt z + , that is, the integrated helicity flux, 
and -£/u nC; box is non-zero. The threshold for Ha/& 2 seems to be —1.7 for w max = 2.0 and —1.4 for 
ttmax = 1-0. w max = 0.5 seems to be too noisy to draw a reliable conclusions; possible reasons for 
that are discussed in the text. 

Fig. [15] compares the generalized twist number to the traditional twist number described 
above. The twist number was computed only for the thinner subvolumes of the torus, 
ttmax = 0.5i? nd zu max = R. Fig. [T5l shows agreement quite well for w max = R and less well 
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for w m ax = 0.5R. The reason might be the following: the smaller the subvolume, the the 
fewer points does it have, so that, first, there are fewer field lines to be traced to measure 
twist, and second, the potential field, obtained by relaxation is numerically less precise. 
Nevertheless, the magnitudes and the general behaviors do agree. 

Fig. [15] also shows the twist number measured for the potential field in a subvolume P, 
is zero to measurement error. Note, that a significant portion of the torus is emerged, its 
length is not large enough (relative to the thickness) for the thin tube approximation to be 
valid. As the twist of the potential field should theoretically be zero (as well as generalized 
twist), this plot also gives an idea of the magnitude of the error of twist measurements; at 
most times the error is less than 15% of the value. 
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Fig. 15. — A comparison between generalized twist number (solid line with diamonds) and the 
"thin tube" classical twist number (dotted line with asterisks) for two subvolumes of a different 
size. Also, the "classical" twist number for a potential field (dashed line with squares). 
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Discussion 



We have demonstrated that, at least in one MHD simulation, the quantity, Tw( gen ), de- 
fined in terms of the additive self helicity shows a threshold beyond which the system became 
dynam ically unstable. The simulation we considered, originally studied by iFan &: Gibson 
( 120031 ). is a three-dimensional, numerical solution of the time- dependent, non-linear evo- 
lution of an emerging flux system. The original study established that the system became 
unstable to a current-driven (kink) mode at some point during its evolution. In this work 
we have shown that the quantity Tw^ gen ) increases until the instability (Tw( gen ) ~ 1.5) at 
which time it drops. This drop occurs as a natural consequence of the instability itself. 

The quantity we propose as havi ng a threshold, Tw( 9Kn ), is c ompu ted using a version 
of the self helicity previous defined by lLongcope fc Malanushenkol (120081 ). The present work 
has provided a detailed method for computing this quantity for any complex bundle of field 
lines within a magnetic field known on a computational grid. We also demonstrate that for 
the very special cases when that bundle can be approximated as a thin flux tube, Tw^ gen ) 
is approximately equal to the traditional twist number, Tw. In the case of thin flux tubes 
which are also dynamically isolated, free magnetic energy is proportional to (Tw) 2 . Their 
free energy may be spontaneously reduced if and when it becomes possible to reduce the 
magnitude of Tw at the expense of the writhe number, Wr, of the tube's axis. 

All this supports the hypothesis that Tw[ gen ) could be treated as a generalization of 
Tw. Such a generalization might be extremely useful in predicting the stability of magnetic 
equilibria sufficiently complex that they cannot be approximated as thin flux tubes. The case 
we studied, of a thick, twisted torus of field lines (IFan &: Gibson! 120031 ). appears to become 
unstable when Tw( gen ) exceeds a threshold value between 1.4 and 1.7. This value happens 
to be similar t o the threshold on T w for uniformly twisted, force-free flux tubes, Tw ~ 1.6, 
as A6 « 3.3tt faood & Priestlll979h . 



Previous investigations have shown that the t hreshold on Tw depe nds on details of the 
equilibrium such as internal current distribution (IHood fc Priestl Il979f ) . It is reasonable to 
expect the same kind of dependance for any threshold on Tw( gen ), so we cannot claim that 
Tw( gen ) < 1.7 for all stable magnetic field configurations. To investigate such a claim is 
probably intractable, but useful insights may be obtained by applying the above analysis to 
magnetic equilibria whose stability to the current-driven instability is already known. The 
paucity of closed-form, three-dimensional equilibria in the literature, and far fewer stability 
analyses of them, suggests this may be a substantial undertaking. 
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